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Abstract 

In biochemical systems some of the chemical species are present with only small numbers of 
molecules. In this situation discrete and stochastic simulation approaches are more relevant 
than continuous and deterministic ones. The fundamental Gillespie's stochastic simulation 
algorithm (SSA) accounts for every reaction event, which occurs with a probability deter- 
mined by the configuration of the system. This approach requires a considerable compu- 
tational effort for models with many reaction channels and chemical species. In order to 
improve efficiency, tau-leaping methods represent multiple firings of each reaction during a 
simulation step by Poisson random variables. For stiff systems the mean of this variable is 
treated implicitly in order to ensure numerical stability. 

This paper develops fully implicit tau-leaping-like algorithms that treat implicitly both 
the mean and the variance of the Poisson variables. The construction is based on adapting 
weakly convergent discretizations of stochastic differential equations to stochastic chemical 
kinetic systems. Theoretical analyses of accuracy and stability of the new methods are 
performed on a standard test problem. Numerical results demonstrate the performance of 
the proposed tau-leaping methods. 
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1. Introduction 

Biological systems are frequently modeled as networks of interacting chemical reactions. 
In systems formed by living cells stochastic effects are very important, as typically some re- 



actions involve only a smal 



Master Equation (CME) 



number of molecules (of one or more species) 



1]. The Chemical 



3| governs the time-evolution of the probability function of the 
system's state. Gillespie proposed the stochastic simulation algorithm (SSA), a Monte Carlo 
approach based on sampling exactly the probability density evolved by the CME 4|. Since 
each reaction is accounted for individually, the overall computational effort becomes an issue 
with systems of practical interest. This motivates the development of approximate sam- 
pling algorithms that trade some accuracy in order to considerably improve computational 
efficiency 

One approximate acceleration procedure is the "tau- leaping method" in which mul- 
tiple reactions are simulated within a pre-selected time interval of length r. The tau-leaping 
method requires that r satisfies the "leap condition" : the expected state change induced by 
the leap must be sufficiently small such that propensity functions remain nearly constant 
during the time step r. In this case the number of times that each reaction fires in the 
interval r is approximated by a Poisson random variable. 

While the tau-leaping method is efficient for single timescale systems, it becomes unstable 
for stiff systems when the stepsize r is large. Stiffness characterizes the dynamics where well- 
separated "fast" and "slow" time scales are present, and the "fast modes" are stable. The 
implicit tau-leaping method improves the numerical stability ^ , but it has a damping effect 
and its results have much smaller variances than SSA results. The trapezoidal tau-leaping 

Q 

formula was proposed to reduce this damping effect [7[. Additional approaches have been 
developed to accelerate the efficiency of the exact SSA through various approximations l], [], 
101 ]. Improved step size (r) selection is discussed in . An alternative point of view is to 



understand the tau-leaping method as the Euler scheme for stochastic differential equations 
(SDEs) 
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131 ]. applied to stochastic chemical kinetics. This is the point of view taken 
in this paper. We propose new tau-leaping-like methods motivated by weakly convergent 
discrete time approximations of stochastic differential equations \\\ . 

The existing implicit tau-leaping methods treat implicitly only the mean part of the 



Poisson variables; the variance part is treated explicitly. Therefore current algorithms can 
be characterized as partially implicit. This paper develops several fully implicit algorithms, 
where both the mean and the variance parts of the random variables are solved implicitly. 
The "BE-BE" method uses the stochastic backward Euler method for both the mean part 
and the variance part of the Poisson variables. The "BE-TR" method uses the implicit 
stochastic trapezoidal method for the variance part of the Poisson variables. The "TR- 
TR" method discretizes both the mean and the variance of the Poisson variables with the 
trapezoidal method. This work also proposes implicit second order weak Taylor tau-leaping 
methods for the stochastic simulation of chemical kinetics. Numerical stability is investigated 
theoretically in the context of the reversible isomerization reaction test problem, an approach 
that is well accepted [3, Q . 

Numerical experiments are performed with three different chemical systems to assess the 
efficiency and accuracy of the new implicit algorithms. The numerical results show that the 
proposed methods are accurate, with an efficiency comparable to that of the original implicit 
tau-leaping methods. They confirm the theoretical stability analysis conclusions that out of 
the six new methods four are unconditionally stable, and two are conditionally stable. These 
analyses perfectly explain our preliminary results reported previously If], [l^ . The numerical 
experiments show that, for stiff systems, all three fully implicit tau-leaping methods avoid 
large damping effects and are stable for any stepsize 16|. But two of the implicit second 
order weak Taylor methods show unstable behavior for large stepsizes (although they are 
more stable than the explicit tau-leaping method 16|). 

The remaining part of the paper is organized as follows. Section[2]describes the traditional 
SSA algorithm. Numerical schemes for the solution of SDEs are presented in Section [3J In 
Section H]the proposed new methods are introduced. Section |5]performs a numerical stability 
analysis using a traditional test example. Results from numerical experiments with three 
different systems are presented in Section |HJ Section [7] draws conclusions and points to future 
work. 
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2. Stochastic Simulation Algorithms for Chemical Kinetics 

In this section we briefly review the traditional SSA and tau- leaping algorithms for 
stochastic chemical kinetics. 

2.1. Exact Stochastic Simulation Algorithm 

Consider a biochemical system involving N molecular species Si, . . ., Sn, composed of M 
reaction channels Ri, . . ., Rm- Denote by Xi(t) the number of molecules of species Si at time 
t. We are interested to generate the evolution of the state vector X[t) = (Xi(t), ...,Xj^(t)) 
starting from an initial state vector X(t ). Assume that the system is well-stirred in a 
constant volume Q and is in thermal equilibrium at some constant temperature. The state 
change vector Vj = u.j = (vij, J^nj) f° r the channel Rj is defined as the change in the 
population of molecule Si caused by one Rj reaction. The propensity function a,j gives the 
probability aj(x)dt that one Rj reaction will occur in the next infinitesimal time interval 
[t,t + dt). 

The SSA simulates every reaction event 4|. With X(t) = x, p(r, j\x,t)dr is defined 
as the probability that the next reaction in the system will occur in the infinitesimal time 
interval [t + r, t + r + dr), and will be an Rj reaction. By letting ao(x) = J2jLi a j( x )> ^ ne 
equation 

p(r,j\x,t) = aj(x)exp(-a (x)T) 

can be obtained. A Monte Carlo method is used to generate r and j. On each step of the 
SSA, two random numbers r\ and r 2 are generated from the uniform (0,1) distribution. From 
probability theory, the time for the next reaction to occur is given by t + r, where 

a (a?) V r i/ 

The next reaction index j is given by the smallest integer satisfying 

j 

y~]aj>(x) > r 2 a (x). 
j>=i 

After r and j are obtained, the system states are updated by X(t + r) := x + Vj, and the 
time is updated by t := t + r. This simulation iteration proceeds until the time t reaches 
the final time. 
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2.2. Tau-Leaping Method 

The SSA is an exact stochastic method for chemical reactions, however, it is very slow 
for many real systems because the SSA simulates only one reaction at one time. One of the 
approximate simulation approach is the tau- leaping method Q. The basic idea of the tau- 
leaping method is that multiple reactions can be simulated at each step with a preselected 
time r. The tau-leaping method requires that the selected r must be small enough to satisfy 
the leap condition, i.e., the expected state change induced by the leap must be sufficiently 
small so that propensity functions remain nearly constant during the time step r. 

Given X(t) = x, denote by Kj(r; x, t) the number of times that reaction channel Rj fires 
during the time interval [t, t + r) where j = 1, . . . , M. The state X(t) = x is updated by 



If the leap condition is satisfied, Kj(r;x,t) can be modeled by a Poisson random variable 
which counts the number of occurrence during a given time period. A Poisson variable 
with parameter a (denoted by V(a)), takes the value k with a probability V(X = k) = 
[e~ a (a) k ]/k\. For stochastic chemical systems V(ar) is interpreted physically as the number 
of events that will occur in any finite time r, given that the probability of an event occurring 
in any future infinitesimal time dt is a dt. Tau-leaping methods use the approximation 



where Vj is a Poisson random variate parameter aj(x)r. 

2.3. Implicit Tau-Leaping and Trapezoidal Methods 

In general, the tau-leaping methods are only able to perform well if they continue to take 
time steps that are of single timescale as fast or slow mode. This drawback is caused by the 
fact that explicit methods advance the solution from one time to the next by approximating 
the slope of the solution curve at or near the beginning of the time interval. For a "stiff" 
system with widely varying dynamic modes among which the fastest mode is stable, the 
leap condition is used to bound the step size r to be within the timescale of the fastest 
mode. Therefore, large leaps are not feasible for stiff systems as they result in no advantage 



M 




(1) 



Kj{r;x,t) « Vj(aj(x)T), 
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compared to the exact SSA. In addition, forced big time step size r might lead to unstable 
population states. 

The tau- leaping method is explicit because the future random state X(t+r) is driven only 
by an explicit function of current state X(t). An implicit tau- leaping method 6| modifies 
the explicit tau-leaping method as follows. Vj can be split as 

We then evaluate the mean value part ajT and the zero-mean random part (variance of the 
Poisson variables) Vj — a^r at the known state X(t). Therefore, 

M 

X (t + r)=x + J2 "i {raj (X(t + r)) + Vj(aj(x)r) - ra,(x)} . (2) 
i=i 

The implicit equation is solved by Newton's iteration method, and the floating point state 
X(t + r) is rounded to the nearest integer values. This implicit tau-leaping method allows 
much bigger step size than the explicit tau-leaping method for stiff systems. But large step 
sizes might provoke damping effect, which means that when a large step size is used to solve 
a stiff system, it yields a much smaller variance and damps out the natural fluctuations of 
the stochastic nature 6|. 

The trapezoidal tau-leaping formula was proposed to reduce the damping effect of the 
implicit tau-leaping formula [7|. The formula is 

M 

X(t + r) = x + {^ a J ( X (t + r )) + ^M^)r) ~ ^(x)} . (3) 

i=i 

Because the trapezoidal rule has a second order convergence without damping effect, this 
formula has better accuracy and stiff stability than the implicit tau-leaping method. The 
trapezoidal method, however, is only second order for the mean value, and still first order 
for the variance. 

3. Discrete Time Approximations for SDEs 

This section discusses the numerical solution of stochastic differential equations (SDEs), 

n 

with an emphasis on weak approximations [14J. 
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3.1. Stochastic Differential Equations (SDEs) 

SDEs are differential equations that incorporate white noise (the "derivative" of a Wiener 
process) and their solutions are random processes. Consider the following d-dimensional SDE 

n 

system [14j 

dX(t) = n(X(t)) dt + a{X(t)) dW(t) , (4) 

X(t) G R d , G R m , t > 0} is an m-dimensional Wiener process, and the functions 

fi : R, d — > R d and a : H d — > H dxm are sufficiently smooth. We call /i the drift coefficient 
and a the diffusion coefficient. 

Because the Wiener process is non-differentiable, special rules of stochastic calculus are 
required when deriving numerical methods for SDEs. There are two widely used versions of 
stochastic calculus, Ito and Stratonovich With Ito calculus, the solution to SDE (TjO) 
can be represented as an Ito integral [l^] 

X(t)=X(t )+ f fi(X(s))ds+ f a(X(s))dW(s), t G [t ,T]. (5) 

With Stratonovich calculus, the solution to (jl]) is 

X(t) = X(t )+ f fi(X{s))ds+ j a(X{s))dW{s), te[t ,T], 

J to J to 

»(X(t)) = ^X{t))-\a{X{t)) d ^-{X{t)), 
where fi is the modified drift coefficient. 

3.2. Convergence 

Consider a time discretization of the SDE ([5]) which uses a maximum step size 5 and 
produces an approximation {Y s (t)} of {X(t)}. The magnitude of the pathwise approximation 
error at a finite terminal time T is measured by the expected absolute value of the difference 
between the Ito process and the approximation jl^] 



e(5) = E[\X(T)-Y\T)\] 



The following two definitions of convergence [14J are useful in the analysis of discretization 
methods. 
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Definition 3.1 (Strong convergence [l^] ) . A time discrete approximation Y 5 {t) with maxi- 
mum step size 5 converges strongly to X at time T if 

\M[\X{T)-Y 5 {T)\] = 0, 

and if there exists a positive constant C , which does not depend on 5, and a finite 5q > 
such that 

E[\X(T) - Y 5 {T)\] < C<F 
for each 5 G (0, 5o), then Y s is said to converge strongly with order 7 > 0. □ 

In many practical situations it is not necessary to have numerical solutions that accurately 
approximate each path of an Ito process. Often one is only interested to accurately compute 
moments, probability densities, or other functionals of the Ito process. The concept of weak 
convergence jl^] describes numerical accuracy in this situation. 

Definition 3.2 (Weak convergence ) . A time discrete approximation Y 5 (t) with maxi- 
mum step size 5 converges weakly to X(t) at time T as 5 j. 0, with respect to a class C of 
polynomials g : M. d — >■ R if 

jm\E[g(X(T))}-E[g(Y s (T))]\=0, 

for all g 6 C. If there exist a positive constant C , which does not depend on 5, and a finite 
5o > such that 

\E[g(X(T))]-E[g(Y\T))]\<C5^ 
for each 5 G (0, So), then Y s is said to converge weakly with order f3 > 0. □ 

These two convergence criteria lead to the development of different discretization schemes. 

3.3. Discretization Schemes 

Consider a time discretization t° < t 1 ■ ■ ■ < t n < ■ ■ ■ < t = T of the time interval [t°, T]. 
The stochastic Euler approximation of the SDE fll]) is 

in 

Y k n+1 = Y k n + fjL k At n + A ^ n , k = 1, • • • , d (6) 

3=1 



8 



where superscripts denote vector and matrix components. We follow our convention in 
writing 

/i fc = /i fc (r,F") and a kJ = a ktj (t n ,Y n ) . 

Here 

AW? = Wf +1 - Wf 

is the N(0; At n ) increment of the jth component of the m-dimensional standard Wiener 
process W on [t n ,t n+1 ], and AW™ and AV^™ are independent for j\ ^ ]2- It was shown [3] 
that the Euler scheme converges with strong order 7 = 0.5 under Lipschitz and bounded 
growth conditions on the coefficients /x and a. 

For weak convergence the random increments AW n of the Wiener process can be replaced 
by other random variables AJW 1 which have similar moment properties to the AH^ n , but 
are less expensive to compute [l^j]. For instance, in the scalar case d = m = 1, a weak Euler 
approximation with weak order (3 = 1.0 is 

y n+i = Y n + fiAt n + a AW n 
where AW™ satisfies moment condition [ijj] 



E 



AW n 



E 



(AW n ) 3 \ I + IE (AW 



rn\2 



At n 



<C(At 



n\2 



(7) 



for some constant C. A simple example of such a random variable is the two-point distributed 
AW n with probability 

P (AW n = ±VAr) = - . (8) 

3.4- The Fully Implicit Euler Scheme 

In the general multi-dimensional case the fcth component of the weak Euler scheme has 
the form 



= n n + At* + 5>*j AW?, Y k ° = X 



(9) 



where AW™ satisfies moment condition (J7J). The family of implicit Euler schemes [14j reads 



yn+l = yn + { a/ife (r+\ Y n+1 ) + (1 - «) ^} At" + ^ <7 fcJ AW? . 



(10) 
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The parameter a here can be interpreted as the degree of implicitness. With a = 1.0 it is 
the implicit Euler scheme, whereas with a = 0.5 it represents a stochastic generalization of 
the trapezoidal method. 

From the definition of Ito stochastic integrals, a meaningful fully implicit Euler scheme 
cannot be constructed by making the diffusion coefficient (a) implicit in an equivalent way to 
the drift coefficient (fi). To obtain a weakly consistent implicit approximation it is necessary 
to appropriately modify the drift term 14j. Such a family of fully implicit stochastic Euler 
schemes is 

Y k n+1 = + {ap k <f+\Y«+ x ) + (1 - a)jll} At" 

m 

+ E {^kA tn+ \ Yn+1 ) + (1 - V)<7kj} A% , (11) 



where AW™ is as in (JS} and the corrected drift coefficient is defined by 

doj_ 

'dx k 



.i 

m d 



j=l k=l 

For a — r\ — 1.0 the scheme (ITT)) is the fully implicit Euler method. For rj = 0.5 the 
corrected drift ffj. = ^ is the corrected drift of the corresponding Stratonovich equation, 
and for a = 0.5 the scheme (TTTj) yields the fully implicit trapezoidal method. 

3.5. The Second Order Weak Taylor Scheme 

In the general multi-dimensional case d,m = 1,2,... the kth component of the second 
order weak Taylor scheme reads j3] 

Y k n+1 = Y k n + A** At" + l -L, » k (At") 2 

m m 

• E{o,, ; Air;' • L a kd I^ ■ /,,//./ '•" } • £ L h a kd2 I^\ (13) 

j = l jlj2 = l 

where operators L and Lj are 



<9t ^ <9x 2 2y f-| ' ' ftr^ dx e ^ ' <9x 2 

2=1 z,£=l h=l z=l 

for j = 1, 2, . . . , m. In addition, the multiple Ito integrals are abbreviated by 



n+l „„2 

e 



10 



Here we have multiple Ito integrals involving different components of the Wiener process, 
which are generally not easy to generate. Therefore f )13p is more of theoretical interest than 
of practical use. However, for weak convergence we can substitute simpler random variables 
for the multiple Ito integrals In this way we obtain from f lT3|) the following simplified 
order two weak Taylor scheme with the kth component 



1 m f 1 1 

y k n+1 = Y k + ^ Ar + -lq ^ (At n ) 2 + + 2 Ar ( L ° a ^ + L J ) 

m 

+ E L n ^ (A% AW? + VU) • 



AH/ n 



(14) 



Jl J2 = l 

Here the for j — 1, 2, . . . , m are independent random variables satisfying moment condi- 
tions 



E[AW r 



E 



(AW 



n\3 



E 



+ 



E 



(AW 



n\2 



At r 



E 



(AW n f 
(AW n f 



?,(At n Y <C(At 



n\3 



(15) 



for some constant C. An N(0; At n ) Gaussian random variable satisfies the moment condition 
(I15p . and so does the three-point distributed AW n with 



1 



P [AW n = ±V3AFJ = -, P (AW n = J 

The Vj lt j 2 are independent two-point distributed random variables with 

1 



for j 2 = 1, • • • ,Ji - 1, 



and 



P(V n , J2 = ±At n ) 



V h , h = -Af 



for j 2 = ji + 1, • • • , m and ji = 1, . . . , m. 



(16) 



(17a) 



(17b) 



(17c) 



4. Implicit Tau-Leaping-Like Schemes 

We now propose several new fully implicit tau-leaping methods motivated by the SDE 
solvers discussed in Sectional 
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4-1. The Fully Implicit Tau- Leaping Methods 

We apply the fully implicit weak Euler scheme (ITT]) to the stochastic chemical kinetic 
problem. Recall the explicit tau- leaping method ((!]). The Poisson variate can be rewritten 
as the mean value part plus the variance part of the Poisson variables. Then the variance 
term is scaled by the standard deviation of dj{x) as below 



Vj(aj(x) t) = a,j(x) r + Jdj(x) AVj 

where the Poisson noise 

Av . = Pj( a j(s) r)-aj(x)r 

3 V a A x ) 

is close to a normal variable iV(0; r) when dj is large. The scheme <^ can be written as 

M M 

X(t + t) = x + ^ Vj aj(x) r + ^uj Ja~(x) AVj . (19) 
i=i j=i 

The weak Euler scheme in vector notation, reads 

m 

yn+1 = Y n + fl At"' + a 3 AW 3 U ( 20 ) 

3=1 

where Oj is the jth column of a. We note that (II 9p is similar to the Euler scheme ( 1201) with 

M 

[i = Vj dj(x) , At n = t , Oj = Vj J dj (x) . (21) 

3=1 

4.1.1. The Fully Implicit "BE-BE" Method 

The fully implicit "BE-BE" tau-leaping method uses the Backward Euler discretization 
for both the mean and variance of the Poisson variables. In (jlip the choice a — rj — 1 
simplifies the fully implicit weak Euler scheme to 

m 

yn+1 = y« + JZ(t n+1 , Y n+1 ) At n + a A tn+1 > Yn+1 ) 

3=1 

where AW™ satisfies moment condition (JTJ). Besides the original random variable AW™ = 
AW™, simpler options like (JS]) are possible fl4|. 

Using (12 ip the corrected drift coefficient ( fT2l can be written as 

M / N 

l ■ 

ji = [i 



l M N 



ddj(x) 
dx k 

3=1 \k=l 
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Finally the "BE-BE" fully implicit tau-leaping method has the form 

m m / N p. 

X(t + r)=x + rJ2 ifli (X(t + r))) - I * £ W + r )) 

3=1 3=1 \k=i k 

M 



+ J2 u i\i a i( X ( t + T )) AW i ( 22 ) 
i=i 

where AW,- = AVj. For large a,-, APj is close to a normal variable and AWj can be replaced 
by a random variable with the correct statistics, e.g., as given by OH]). 

4.1.2. The Fully Implicit "TR-TR" Method 

The fully implicit "TR-TR" method uses an implicit trapezoidal discretization for both 
the mean of and the variance of the Poisson variables. The choice a = rj = 0.5 in ffTTT) leads 
to 

Y n+1 = Y n + - {]J(t n+ \ Y n+1 ) + 71} At n + - EM r+1 ' yn+1 ) + a A AW i » 

3=1 

where the corrected drift coefficient (fT2|) is 

., m d 

^ = ^2LL^' (23) 

3=1 k=l 

and is equivalent to the Stratonovich drift coefficient /x. 

From (12 ip the "TR-TR" fully implicit tau-leaping method has the form 

M 



X(t + r) = x+ T -Y, "j («i (X(t + r)) + a,(x)) 



3=1 

M f 1 w 

-^E^izE 



dcij(X(t + r)) da,j(x) 



j=i I fc=i re 



+ 2 £^ (v / ^( x ( t+r )) + v^)) 
3=1 ^ ^ 

where the AWj = AT^- or, for large a,, can be replaced by ([8]). 



(24) 



The Fully Implicit "BE-TR" Method 
The fully implicit "BE-TR" method uses a backward Euler discretization for the mean 
(deterministic) part, and the implicit trapezoidal discretization for the variance. In (fTTl) the 
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choice a = 1.0 and 77 = 0.5 simplifies the fully implicit weak Euler scheme to 

Y n+i =Y n + /j (t n+ \ Y n+1 ) At n + - ^{0i(t n+ \ r ' i+1 ) + (T i(^ ) ^n)} AWj , 

i=i 

where the corrected drift coefficient ([12]) is equal to ([23]). From Q2D the "BE-TR" fully 
implicit tau-leaping method has the form 

X{t + r) = x + t £ „, + r)) - I £ „, f £ ^ ftbW + r)) 



3=1 j=i \k=i 



+ g E "i ( y/^(X(t + T )) + \fe ) A ^ (25) 
j= i V / 

where the AWj = AVj or, for large be replaced by (jSJ). 

Implicit Second Order Weak Taylor Tau-Leaping Methods 

The simplified order two weak Taylor scheme (JHJ) motivates the following family of 
methods for stochastic kinetic equations: 

Y k n+1 = Y k n + {aft k (t n+1 , Y n+1 ) + (1 - a) fi k } At n 

+ 1(1 - 2a) {(3L of i k (t n+1 ,Y n+1 ) + (1 - /3)L /i fe } (At n ) 2 



Jl = lj2 = l 

m f 1 1 

+ E a ^ + 9 (L ° °^ + (1 _ 2a)L ^ ^ Ar r • 
.7=1 ^ ^ 



(26) 



4-2.1. Implicit Second Order Weak SSA with a = 1.0 and (3 = 1.0 
When a = 1.0 and /3 = 1.0 the scheme f[26]) becomes 



Y k n+1 = Y k n + fi k {t n+1 ,Y n+1 )At n - h 0f i k {t n+1 ,Y n+1 ){At n ) 2 



■ \ E ( A »7 AlT X • V h,i 

ii=lj2=l 

m r 1 1 

+ E [ ff *j + 2 (Lo ^ " L ' ^ fc) Ar / Aw ^ n • (27) 
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We apply the implicit order two weak Taylor scheme to the stochastic chemical kinetic 
problem in a similar manner to the fully implicit tau-leaping methods. Note that 

k=l k,l=l h=l K * k=l 

d ri -, d m O o d ri 

v-^ e'er,- 1 CTj , r V"^ (XT,-, , 00 , 

= L%+2LL^ ff ^ and ^ = ^^—.(28) 

fc=l K fc/=l h=l K £ fe=l ft 

From (I2T1) . (j27|) . and ( 128]) the implicit order two weak tau-leaping SSA method with a = 1.0 
and /? = 1.0 has the form 

M 

X{t + r) = x + tJ2 v i ( a j ( x (t + r ))) 



yE^ E — ^ — (2>w 

j=i k fc=i \/i=i 



1 ^E d 2 a,(X(t + r)) . 



+ 1 E ^™ I E ( E ^ + 



M ( N / M ~ I \ 

+E ^\A^)-^\A^)E^- E^l^lh^ 

i=i I fc=i \/i=i fc 

M ( N / s / M 

T \ J \ ^ oaj{x) 



31,32 



*jw M=1 — « \ h=1 

4-2.2. Implicit Second Order Weak SSA with a = 1.0 and /3 = 0.0 
When a = 1.0 and /3 = 0.0 the scheme (1261) reads 

n n+1 = n n + Vk(t n+ \ Y n+1 )At n - h /i fc (Ar) 2 



(29) 



+ 1 E /-,: n;,, (aiT;;' Ail;;; • \ },.,-/ 

jl=lj2=l 
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The corresponding implicit order two weak tau-leaping SSA method has the form 

M 

X {t + r)=x + rY, Vj (X(t + r))) 
3=1 

M ( N n , s / M \ „ N no / \ / M 



j=i Lfe=i K \h=i / fc,<=i fc * \h=i 

+ \ E ^7^75 ( E \/^> f E <^^r) ( aR7 » a ^ + *5.. 



j=l k k=l \h=l 

M ( N a , , / M 

t v— «, Vj I r-v oaj{x) 



.i 



4 U- W a i( x ) I fc =i \ h =i 

d 2 aAx) ' 



^P^^l' 4 (30) 



4-2.3. Implicit Second Order Weak SSA with a = 0.5 

When a = 0.5 the scheme (126]) does not depend on /3. The method reads 



+ 2 E rr,,,, (aiT;;; Ail;;; • \} : . ; ,) 

jl=lj2 = l 



J-- 

The implicit order two weak tau-leaping SSA method for a = 0.5 has the form 

M M 

X{t + r) = x+ T -Y^ uj {a, (X(t + r)) + a,(x)} + ^ Uj \Ja~{x) AWj 

3=1 3=1 



M ( N a , s. / M 

4^ 4v ^M\^ dx k \fe ^ 



(31) 
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5. Stability Analysis 



In this section we perform a theoretical stability analysis of the fully implicit methods 
proposed in Section HI Specifically, we take the well established approach [l5|, llO] of applying 
the methods to the reversible isomerization model and comparing the discrete results with 
the available analytical solution. 

5.1. Reversible Isomerization Model 

Following Rathinam et al, flol [ic| we consider the reversible isomerization reaction 



system 



Si i S 2 . 

C2 



(32) 



Let X t denote the population (number of molecules) of Si at time t, X T the total population 
of 5i and 62, and 

A = ci + c 2 . (33) 

Usually the case with ci = c 2 is considered. Note that X T is constant in time, and therefore 
the population of S 2 at time t is X T — X t . The deterministic reaction rate equation for this 
system is the ODE: 

dX t 



dt 



ciX t + c 2 (X 7 - X t ) = -XX t + C2X 1 . 



Therefore the mean E,[X t ] and variance Var[X t ] satisfy the following ODEs: 

<mx t ]__ XE[Xt]+C2XT 



dt 

dVax[X t ] 
~dt 



-2A VarLY f ] + c 2 X T + ( Cl - c 2 ) E[X t 



As t goes to infinity, the asymptotic value of the exact mean E[X^] and the exact variance 
Var[X^] are [l5|,[l3] 

(34) 



TO = VarK] = ' 



A ' -- l -°° J A 2 ' 
5.2. Stability Analysis of the Traditional Tau-leaping Methods 

Recall the explicit tau-leaping method ([T|). Applying the explicit tau-leaping method 
with a fixed step size r to the test problem (132]) gives 

X n+ i =X n - Vi(cirX n ) + V 2 (c 2 r(X T - X n )) , (35) 
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where X n is the numerical approximation of X t at time t n . 



The following lemma about the conditional probability from 19| will prove useful for the 
derivation. 

Lemma 5.1. If X and Y are random variables, then 

E[Y] = E[E[Y\X}}, 
Var[F] = E[Var[F | X]] + Var[E[F | X]]. 

By Lemma I5.1[ the mean of the Eq. (|35|) is 

E[X n+1 ] = (1 - Ar) E[X n ] + c 2 X T r . 

This imposes the stability condition 

|1-At|<1, (36) 

which implies < Ar < 2 for the stepsize. For n — > oo we obtain the asymptotic mean 

c 2 X T 



E[X C 



A 



For the variance we have 



Var[X n+1 ] = (1 - Ar) 2 Var[X n ] + ( Cl - c 2 ) rE[X n ] + c 2 X T r . 



(37) 



The stable domain for the variance is given by (1 — Ar) < 1 and is the same as ( I36p . For 
n — >■ oo in (I3T|) . the asymptotic variance is 

2 



Var[X c 



2 - Ar 

Thus the variance given by the explicit tau-leaping method does not converge to the theo- 
retical value, even if the stability condition is satisfied. If Eq. ( 136]) is satisfied, VarLXoo] is 
larger than VarLY^]. 

Similarly, the stability region, asymptotic mean, and asymptotic variance for the tradi- 
tional implicit tau-leaping method are 



< 1, E[X, 



1 + Ar 

For the trapezoidal method 
2 - Ar 



oo^^ = E[X*J, Var[X 



2 + A^ 



VarK 



2 + Ar 



<1, E[X OQ ] = ^=E[X* OQ ], VarlXoo] = = Var[X* 



A 



A 2 



(3J 



(39) 
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5. 3. Stability Analysis of the Fully Implicit Tau-Leaping Methods 
Recall the BE-BE fully implicit formula 



X(t + r) = x + £ + r)) - \ (± ^ ^ ( ^ +T) ) 



+ J aj (X(t + T))T 



Vj(a,j(x) r) — a>j(x) t 



V a j( x ) 

We apply the BE-BE tau-leaping methods with a fixed step size r to the test problem 
For N — 1, M = 2, = —1, z/i i2 = 1, ai(a;) = ciX, and 02(0;) = c 2 (X T — X), we have that 



ra+1 



X n - rAX„ +1 + r ( c 2 X T - ^ + % 



2 2 
Pi(rciX n ) - rciX, 



, / YT y ( ^(rc 2 (X^-X n ))-rc 2 (X y -X w ) \ 
v n+1 l V* r -X n J 



(40a) 
(40b) 

(40c) 



Derivation of the mean for the simplified equation ( 1401) is quite intricate due to the square 
root in the denominator. In order to derive the stability region we first employ an inequality 
condition. Denote by E n [-] = E[-|X„]; from lemma O E[ ■ ] = E[E n [ • ]]. Taking the 
expectation of f!40bp leads to 



E K 



X n +\ 



VxircxXn) - TdX n 



< - E n [X n+1 ] + X - E n 



(Vl(T Cl X n ) -TCiX n ) 

X„ 



E n [X n 



+1 



;i~ci, 



which implies that 



E K 



+ 1 



Pi(rciX n ) - rciX n 



< -E[X n+1 ] + -r Cl . 



(41a) 



Similarly, the expectation of (j40cj) satisfies 

V 2 (rc 2 (X T - X n )) - rc 2 (X T - X n ) 



E 



+ 1 



y/x* - x n 



< E[X T -X n+1 ]+-TC 2 . 



(41b) 
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Plugging (14 lap and (]41bj) into (140]) and taking E[ • ] gives 



E[X n+1 ] < E[X n ] - r\E[X n+1 ] + r (c 2 X T - | + | ) 



2 2 

i E [X n+1 ] + i rci + J E [X T - X n+1 ] + i rc 2 , 



which can be simplified to 



E[X 



n+l. 



< 



(1 + Ar) 

This imposes the sufficient stability condition 

1 



2rc 2 + 2tc 2 X t + X T 
(2 + 2Ar) ' 



< 1. 



(42) 



(43) 



1 + At 

The second approach for the stability analysis is using the Poisson approximation method. 
Recall that the Poisson random variable can be rewritten as the mean value plus the random 
deviation from the mean part 



Vj(aj(x) r) = aj {x)T + yjaj{x) AVj. 



If Oj is large the Poisson noise AVj is close to a normal variable N(0; r). In this case the 
Poisson variable with mean a,j(X(t + r)) r can be approximated by 



V( aj (X(t + r)) r) « aj (X(t + t))t + Ja,(X(t + r)) AVj . 



(44) 



With this approximation the "BE-BE" fully implicit method has the alternative form 

X(t + r) = x + Y t ^MX(t + r)) r) - \ g ^ ^ • (45) 

Applying the alternative BE-BE formula ( 145]) with a fixed step size r to the test problem (132]) 
gives 



X n+1 = X n - Pi( Cl rX n+1 ) + P 2 (c 2 r(X T - X n+1 )) - ^( Cl - c 2 ) 
Denoting by E n+ i[ • ] = E[ ■ |X n+ i] and taking E n+ i of ( 145]) leads to 

X n+l = E n+1 [X n ] - cirX n+1 + c 2 t(X t - X n+l ) - -(ci - c 2 ), 



(46) 



i.e., 



E n+ i[X n ] = (1 + Xr)X n+1 - c 2 rX L + -(ci - c 2 



(47) 
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Then by Lemma 15.11 we have 



E[X n ] = E[E n+1 [X n ] } = (1 + Ar) E[X n+1 ] - c 2 rX T + -( Cl - c 2 ). 



Therefore 



which imposes the stability condition 



1 



< I. 



(48) 



(49) 



1 + Ar 

This approximate stability region is same to the sufficient BE-BE stability condition (1431) 
calculated via inequalities. We conclude that the BE-BE stability is similar to that of the 
traditional implicit tau-leaping method for the reversible isomerization test model. 

The Poisson approximation (|44p allows to deduce the asymptotic mean and variance of 
the approximate solutions 045p . Letting n — > oo in fT45|) we obtain 

nx^] = I (c 2 x t + 

For ci = C2 (the common setting of the test problem) 

E[X O0 } = ^=E[X* oo }. 
The conditional variance of (146]) with respect to X n+ i is 



Var[X n |X n+ i] = (c 2 - Ci)rX n+l - c 2 tX . 



Therefore 

E[Var[X n |X n+1 ] ] = (c 2 - Cl )rE[X n+1 ] - c 2 tX t . 
The variance of ( )47|) is 



(50) 



Var[E[X n |X, 



n+lj 



;i + Ar) 2 Var[X n+1 ]. 



(51) 



From Lemma 15.11 ( )50|) , and ([51 



Var[X n ] = (1 + Ar) 2 Var[X n+1 ] + (c 2 - Cl )r E[X n+1 ] - c 2 rX q 
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Method 



BE-BE 
TR-TR 
BE-TR 



Stability condition 



1- 


fAr 


2- 


-At 


2- 


rAr 




1 


1- 


fAr 



< 1 

< 1 

< 1 



E[X C 



Var[X c 



5?Var[Xi 



2+At 

Var[X^] 



Table 1: Behavior of fully implicit methods applied to the reversible isomcrization problem. All methods are 
unconditionally stable and yield the exact asymptotic mean. TR-TR provides the exact asymptotic variance 
as well. 

Letting n — >■ oo 

VarLYj = (1 + Ar) 2 Var[X 00 ] + (c 2 - c^rE^] - c 2 rX T . 
After replacing the E[X C 



4Ux t ■ Ci ~ C2 



For ci = c 2 as the E[X C 
VarLYj 



A 

Var[X Q 
2 Cl c 2 X T 



Acic 2 X T + (ci - c 2 ) 
2A 2 (2 + Ar) 



CiC 2 X^ 



Var[X^]. 



A 2 (2 + Ar) 2 + Ar A 2 2 + Ar 

This asymptotic variance of the approximate BE-BE (|22|) is same as that of the traditional 
implicit tau-leaping method (138]) . 

A similar approach can be used to obtain the stability region, the asymptotic mean, 
and the asymptotic variance of the TR-TR ( liHj) and BE-TR fl25|) methods. The results are 
summarized in Table [TJ 

5.4- Stability Analysis of the Implicit Second Order Tau-Leaping Methods 

Application of the implicit second order method with a = 1.0 and (3 = 1.0 ( 129]) to the 
test problem (152]) yields 

X n+1 = X n + t(c 2 X t - XX n+1 ) + ~ (n - r 2 - r 3 + r 4 ) + r 5 + r 6 + ^-(c 2 ^ T - AX n ), (52) 
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with 

{Pi(rciX n ) - rciX n } 2 

n = tp h ciV%i, 

{P 2 (rc 2 (X T - X n )) - tc 2 {X t - X n )} 2 
ri = X T -X + ' 2 ' 



r 3 



{^(tcM - r Cl X n } ■ {V 2 {tc 2 {X t - X n )) - rc 2 (X T - X n )} Cl c 2 (X T - X n ) 



Xn 



X n 



V 2 A, 



{Vi(rc x X n ) - r Cl X n } ■ {V 2 (rc 2 (X T - X n )) - rc 2 (X T - X n )} [^X, 



X T - Xn 



X T -X 



r 5 =(l+ y) {V 2 (tc 2 (X t - X n )) - tc 2 (X t - X n ) - V 1 (rc 1 X n ) + r Cl X n ) 
(XX n - c 2 X 



.ViiraXJ-raXn , V 2 (tc 2 (X t - X n )) - rc 2 (X T - X n 



X n X T — X n 

where The Vj lt j 3 are independent two-point distributed random variables as (1171) . In order 
to derive the mean of equation (1521 . we first compute E n [ri], ...,E n [r 6 ]. Using E n [Vi 5 i] = — t, 



En[n] =E, 



X n 



ciVl, 



Var (V^raXn)) 

X n 



TCi = . 



Similarly, E n [r_j] = for j = 2, . . . , 6. Therefore 



(l + Ar)E n [X n+ i] = 1 



AT \ I \t 

— \E n [X n ]+Tc 2 X T ll + — 



From Lemma 15.11 the mean of the numerical solution satisfies 

E[X„ +1 ] = 



2 + 2Ar 

which implies the stability restriction 

2 - A 2 t 2 



2 + 2\t 



(53) 



< 1 0<Ar<l + v / 5. 



(54) 



2 + 2Xt 

The second order weak Taylor method with a = 1.0 and /3 = 1.0 is conditionally stable. For 
the asymptotic mean of the second order weak Taylor method with a = 1.0 and /3 = 1.0, let 
n — > oo in fl53|) . Then we obtain 

c 2 X T 



E[X 



A 



(55) 



which is equal to its exact value 
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The stability condition and the asymptotic mean for the implicit second order with 
a = 1.0 and (3 = 0.0 (I5U|) are calculated in a similar manner, and the results are the same 
as ([MD and fl55|) . 

Application of the implicit second order method with a = 0.5 ( l3Tj) to the test problem fl32l) 
gives 

X n+1 = X n + ^{2c 2 X T - XX n+1 - \X n ) + - (n - r 2 - r 3 + r 4 ) + r 5 + r 6 , (56) 

with 

{V^Xn) - r Cl X n } 2 

n = T? 1" ClVi,!, 

{^ 2 (rc 2 (X T - X„)) - rc 2 (X T - X n )} 2 

r 2 = ~ TPf. 7? — + C 2 V2,2, 



r 3 



r 4 



gi(rcA) - r Cl X n } ■ {V 2 (tc 2 (X t - X n )) - rc 2 (X T - X n )} 

X„ 



l c 1 c 2 {X T - X n 

x„ 



-V 2 ,i, 



{Pi(r Cl X n ) - r Cl Xn} ■ {V 2 (rc 2 (X T - X n )) - rc 2 (X T - X n )} I Cl c 2 X n 



x T -x n 

r 5 = V 2 (tc 2 (X t - X n )) - rc 2 (X T - X n ) - V 1 (rc 1 X n ) + TCl X n , 



X T -X 



■Vi, 



r6= 16 



(XX n - c 2 X' 



T .( V^tcM - TCl X n t V 2 {rc 2 {X T -X n ))-rc 2 {X T -X n 



X„ 



x T -x„ 



Similar to the calculation for the implicit second order weak SSA with a = 1.0 and /3 = 1.0, 
taking expected value E n and then E gives 



2 - \i 



2 + Xt 

The asymptotic stability of E[X n ] requires 

2 - At 

< 1 



E[X„] 



2rc 2 X r 
2 + At 



(57) 



< Ar. 



(51 



2 + Ar 

Because Ar is always greater than zero, the second order weak Taylor methods with a = 0.5 
is unconditionally stable. The condition (158]) is the same as that ( )39|) of the trapezoidal 
tau-leaping method. Letting n — > 00 we have 

c 2 X T 



E[X C 



A 
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which is equal to its exact value (p?j 

Deriving analytically the asymptotic variances for the second order weak Taylor methods 
becomes a very intricate task. For the variance of the implicit second order method with 
a = 0.5 fl3Tj) to the test problem fl32|) . we still use the fact 

Var[X n+1 ] = E[ Var[X n+1 |X n ] ] + Var[E[X n+1 |X n ] ] 

using Lemma (15. II) . By (j57j) . 

2-Ar x 2 



Var[E[X n+1 |X n ]] = ( 1 Var[X n ] 

To calculate the term E[ Var[X n+1 |X n ] ], we should consider the expectation of the variance 
of ( 1561) . This involves the estimation of E[^S-] and E[ xT j_ x ] which cannot be obtained 
simply. This intractable calculation will be analyzed in future work. 

6. Experimental Results 

This section presents numerical results for the new implicit tau-leaping methods applied 
to three different systems. A fixed stepsize strategy is used in each simulation for all methods; 
this allows for a clean comparison of the performance of different algorithms. 

6.1. The Decaying-Dimerizing Reaction Set 

The decaying-dimerizing system 
reactions 

S 1 0, 

S 1 + S 1 ^ S 2 , (59) 

C3 

Q c 4 v Q 
D2 > 03. 

We chose the following values for the parameters 

Cl = 1, c 2 = 10, c 3 = 1000, c 4 = 0.1, 



10j consists of three species Si, S2, and 53 and four 



which will render the problem stiff. The propensity functions are 

ai = X x , a 2 = 5Xx(Xx - 1), a 3 = 1000X 2 , a 4 = 0.1X 2 , 
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Dimer species trajectories by the SSA 
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Figure 1: Time evolution of the numbers of molecules in the decaying-dimerizing problem (|59|). The simu- 
lation is carried out using Gillespie's SSA method. 

where Xi denotes the number of molecules of species Si. The initial conditions are 

Xi(0) = 400, X 2 (0) = 798, X 3 (0) = [molecules]. 

The final time is T = 0.2 seconds. Figure Q] shows the species evolution for the reaction 
set fl59l) solved with the original SSA. 

In order to compare the solutions given by different methods we consider histograms 
of Xi, the number of molecules of Si, at the final time T = 0.2 seconds. Specifically, an 
ensemble of simulation results is carried out for each method, and the final distribution of 
the numerical X\ is plotted as a histogram from 100,000 independent simulations. 

Figure [2(a) shows the histograms of Xi for the decaying-dimerizing system (1591) simulated 
with Gillespie's SSA and with the traditional explicit tau-leaping, implicit tau-leaping, and 
trapezoidal tau-leaping methods. A fixed stepsize r = 2 x 10 -4 seconds is used. Figure 12(b) 
also shows the histograms generated with Gillespie's SSA, and with the methods proposed 
herein: fully implicit BE-BE, TR-TR, BE-TR, implicit order two weak Taylor with a = 1.0 



and (3 = 1.0, a = 1.0 and (3 = 0.0, and a = 0.5. The same fixed stepsize r = 2 x 10" 
used. 



is 



Figures [21 (a) and (b) reveal that the histograms of the trapezoidal tau-leaping method, 
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(a) Histograms obtained with Gillespie's SSA, and (b) Histograms obtained with Gillespie's SSA, the 
with the traditional explicit, implicit, and trapezoidal new fully implicit methods, and the new implicit order 
tau-leaping methods. two weak Taylor tau-leaping methods. 

Figure 2: The histograms of the number of molecules X\ at the final time for the decaying-dimerizing reaction 
system ([59]) . All histograms are based on 100,000 runs of the corresponding methods with a fixed stepsizc 
t = 2 x 10~ 4 seconds. 



fully implicit TR-TR method, and implicit order two weak Taylor method with a = 0.5 are 
closer to the reference (SSA) histogram than those of other methods, for the specific time 
step chosen. 

The explicit method gives very unstable and varying results. Other implicit order two 
weak Taylor methods with a = 1.0 provoke a little wide varying results, but those escape 
the damping effect such as implicit tau-leaping method in Figure |2] (a). From the stability 
analysis, we have proved that the implicit order two weak Taylor methods with a = 1.0 are 
unstable for large stepsizes, and these experimental results confirm the conditional stability. 

In order to numerically assess the accuracy of each method, we carry out simulations with 
different stepsizes, and obtain the corresponding histograms. For each method and step size 
the numerical errors are quantified by the difference between the numerical histograms and 
the reference (SSA) histogram. Two metrics of the difference are employed: the Kullback- 
Leibler (K-L) divergence 20] and the distance metric. 
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The K-L divergence is a non-commutative measure of the difference between two probabil- 
ity distributions P and Q, typically P representing the "true" distribution and Q representing 
arbitrary probability distribution. Therefore we set P to be the distribution obtained from 
SSA, and Q the distribution obtained with one of the other formulae. The K-L divergence 
is defined to be 



where Q(i) 7^ 0, and the summation is taken over the histogram bins. Smaller values of 
K-L divergence represent more similar distributions. Because K-L divergence is not useful 
when there exists zeros for Q, we also use the distance metric, which measures the difference 
between two distributions by 



Here AX is the bin size of the histogram. 

Table |2] shows these metrics based on 100,000 samples generated by different methods 
for fixed stepsizes r = (8/k) x 10 -4 where k = 1,2,4,8. The results show that the mean 
is accurately computed by all accelerated methods. However, the variance and distance are 
different for each formula. For example, the explicit tau formula becomes very unstable for 
a stepsize of 4 x 10 -4 seconds. The implicit tau-leaping, BE-BE, BE TR are far superior to 
explicit tau, but those formulae produce smaller variances compared to the variance of the 
exact SSA that is called as damping effect. 

Three methods (the trapezoidal-tau, the fully implicit TR-TR, and the implicit second 
order weak Taylor with a = 0.5) generate accurate variance results even with large stepsizes. 
The fully implicit TR-TR results are the most accurate among all methods for similar time 
steps, as demonstrated by the smaller distance to the reference histogram in Table [2j The 
implicit second order weak Taylor methods with a = 1.0 are accurate until they become 
unstable for large stepsizes. 

The elapsed CPU times for each method are presented in Table El Figure [3] considers the 
relationship between accuracy and computation time for each of the accelerated methods. 
From the figure, the trapezoidal tau-leaping, the fully implicit TR-TR, and the implicit 
second order weak Taylor with a = 0.5 methods generate accurate solutions with a large step 
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Tabic 2: The mean, variance, K-L divergence, and distance for X\ at T = 0.2 based on 100,000 samples for 
different stcpsizes of the decaying-dimerizing reaction system (f59l) . 
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oo 


oo 


0.080 
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oo 


CO 
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Mean 
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weak Taylor 


Variance 


356.93 


350.17 


348.72 


348.89 


(a = 0.5) 


K-L div. 


0.004 


0.003 


0.002 


0.002 




Distance 


0.625 


0.421 


0.386 


0.318 



size (r = 8 x 10~ 4 seconds) and in a short CPU time. For comparison, 100,000 simulations 
using the SSA took 16,210 CPU seconds, while 100,000 simulations of the fully implicit 
TR-TR took only 377 seconds (2.3% of the SSA time) and provided an accurate solution 
(distance value is only 0.276). The implicit second order weak Taylor method of the a = 0.5 
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Table 3: Elapsed CPU times (in seconds) for each method and time step for 100,000 simulations of the 
decaying-dimerizing reaction system (1591) . 



CPU time (seconds) 


Stepsize (r in seconds) 


Method 


8 x 10~ 4 


4 x 10~ 4 2 x 10~ 4 1 x 10~ 4 


Gillespie SSA 


16210.13 


Explicit tau-leaping 


27.32 


46.91 


130.55 


260.24 


Implicit tau-leaping 


170.57 


340.58 


657.51 


1389.29 


Trapezoidal tau-leaping 


180.42 


350.66 


688.98 


1301.21 


Fully implicit BE-BE 


344.98 


686.49 


1395.1 


2638.74 


Fully implicit TR-TR 


377.06 


746.24 


1400.96 


2752.39 


Fully implicit BE-TR 


340.65 


690.56 


1373.31 


2657.25 


Implicit 2.0 weak Taylor (a = l,/3 = 1) 


398.23 


784.43 


1587.69 


3121.32 


Implicit 2.0 weak Taylor (a = 1, @ = 0) 


391.31 


765.39 


1532.98 


3076.23 


Implicit 2.0 weak Taylor (a = 0.5) 


381.34 


752.84 


1425.83 


2798.54 
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Figure 3: Relationship between solution accuracy (measured by the distance (|6Tj) between the accelerated 
method and the SSA produced histograms) and CPU time for different methods applied to the decaying- 
dimerizing reaction system ([591) . 



with r = 8 x 10 4 fixed step took 381 seconds and produced results of similar accuracy. 



6.2. Schldgl Reaction Set 



The Schldgl reaction model 15[ is a simple but famous bistable system. The system 
contains four reactions 



(62) 



B\ + 2S ¥^ 3S, 
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(a) Histograms obtained with Gillespie's SSA, and (b) Histograms obtained with Gillespie's SSA, the 
with the traditional explicit, implicit, and trapezoidal new fully implicit methods, and the new implicit order 
tau-leaping methods. two weak Taylor tau-leaping methods. 



Figure 4: The histograms of the number of molecules X at the final time for the Schlogl bistable system (|62[). 
All histograms arc based on 100,000 runs of the corresponding methods with a fixed stepsize r = 0.4 seconds. 



where B 1 and B 2 are buffered species whose populations are assumed to remain constant 
over the time interval. 

a = 3 x lcr 7 , c 2 = lcr 4 , c 3 = icr 3 , c 4 = 3.5, JV 1 = ix io 5 , iv 2 = 2x io 5 . 



which will render the bistable system. Hence the propensity functions are given by 

a 1 = ^N 1 X(X-l), a 2 = ^X(X-l)(X-2), a 3 = c 3 iV 2 , a 4 = c 4 X 
2 o 

where X denotes the number of molecules of species S. Initial condition X(0) = 250 at 
T = 0, and final time T = 4 second. 

The histograms generated from 100,000 independent samples of SSA, existing improved 
SSA methods, and proposed methods including fully implicit tau-leaping methods and im- 
plicit order two weak Taylor methods with fixed stepsize r = 0.4 are shown in Figure HI 
We notice that the histogram given by the trapezoidal tau-leaping method, fully implicit 
TR-TR method, and implicit order two weak Taylor method with a = 0.5 are very close 
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Table 4: The mean, variance, distance, and elapsed CPU times (in seconds) for X at T = 4 based on 100,000 
samples for different stepsizes of the Schlogl bistable system (|62|) . 





Stepsize (r in seconds) 


Method 


Metrics 


0.8 


0.4 


0.2 


0.1 


Gillespie 
SSA 


Mean [ Var J 
CPU time 


305.2 (46465.9) 
682.96 


Explicit 

ut*LL — ItJcllJlIlg 


Mean (Var) 

Til Q"f" CI Yl OC* 

CPU time 


296.9 (40957.6) 
5.680 
1.41 


306.2 (42915.6) 
3.155 
2.1 


309.5 (44981.5) 
2.057 
3.43 


308.5 (45929.9) 
1.860 
6.21 


Implicit 

bdlX ~ 1CCIU1I1K 


Mean (Var) 

Tii Q-f- anpn 
LdillLL. 

CPU time 


343.4 (52245.0) 

A A(ZA 

4.404 

4.41 


326.3 (49876.8) 

O Q77 

z.o ( ( 

7.03 


316.9 (48364.8) 

O 1 

z.loD 

12.24 


315.1 (47644.7) 
i.yoD 
22.4 


Trapezoidal 
taudeaping 


Mean (Var) 
Distance 
CPU time 


324.6 (47837.6) 
z.uou 
4.2 


317.4 (47161.6) 
6.79 


312.6 (46727.0) 

X .O^ttJ 

12.07 


311.2 (46719.1) 

X .oxo 

22.6 


Fully implicit 

1)1. 1)1. 


Mean (Var) 

Til Qf" QT1PP 

CPU time 


316.4 (51137.7) 
4.360 
8.64 


318.8 (49359.6) 
2.808 
13.6 


313.5 (47919.2) 
2.158 
23.74 


312.2 (47401.1) 
1.956 
43.86 


Fully implicit 

X XX J. XV 


Mean (Var) 

1 lief onpn 

CPU time 


316.2 (47195.7) 
1.943 
8.13 


312.4 (46743.9) 
1.857 
13.63 


312.2 (46624.0) 
1.836 
24.51 


309.9 (46601.9) 
1.818 
46.4 


Fully implicit 
BE-TR 


Mean (Var) 
CPU time 


335.5 (51920.4) 
4.417 
8.80 


322.3 (49566.8) 
2.761 
13.38 


315.9 (48011.9) 
2.147 
24.98 


311.1 (47325.1) 
1.917 
46.76 


Implicit 2.0 
weak Taylor 
(a = 1,0 = 1) 


Mean (Var) 
Distance 
CPU time 


1122.4 (51112.5) 
3.501 
12.53 


310.3 (49157.7) 
1.890 
18.72 


310.2 (47332.8) 
1.830 
30.98 


310.0 (46612.9) 
1.766 
55.08 


Implicit 2.0 
weak Taylor 
(a = 1,0 = 0) 


Mean (Var) 
Distance 
CPU time 


296.4 (50810.1) 
2.475 
11.74 


306.2 (46870.6) 
1.869 
17.48 


309.5 (46566.0) 
1.842 
28.76 


309.7 (46498.5) 
1.839 
52.64 


Implicit 2.0 
weak Taylor 
(a = 0.5) 


Mean (Var) 
Distance 
CPU time 


313.2 (47441.4) 
1.862 
10.71 


309.9 (46880.3) 
1.840 
16.34 


309.7 (46494.3) 
1.809 
26.47 


310.2 (46503.7) 
1.803 
50.23 



to the exact SSA method than other methods for the specific time step as the histogram 
of the decaying-dimerizing system. The histograms produced by the fully implicit BE-BE 
and BE-TR exhibit damping effect (sharp peaks) while the histogram given by the implicit 
order two weak Taylor method with a = 1.0, /3 = 1.0 and a = 1.0, = 0.0 methods provoke 
a little wide varying results (broad peaks). 

Table H] shows the mean, variance, distance, and elapsed CPU times based on 100,000 
samples generated by different methods for fixed stepsizes. Four fixed stepsizes r = 0.8 /k 
where k = 1,2,4,8 were selected to evaluate accuracy for each time step. The variance 
for all methods are large for the bistability property of the system. Proposed fully implicit 
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TR-TR, and the implicit second order weak Taylor with a = 0.5 produce accurate results 
even with large stepsize r = 0.8. 
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Figure 5: Relationship between solution accuracy measured by the distribution distance (|61j) and CPU time 
for different methods applied to the Schlogl bistable system ([52)) . 



Figure [5] shows the relationship between distance of two distributions (the SSA and each 
accelerated method distributions) and computation time for the different stepsizes of Schlogl 
bistable system. As the previous dimer reaction system, the fully implicit TR-TR and the 
implicit second order weak Taylor method with the a = 0.5 show small distance (good 
accuracy) compared to other accelerated methods with the big stepsize r = 0.8. 100,000 
simulations of the fully implicit TR-TR method with the r = 0.8 took 8.13 seconds with 
accuracy. With the limited results investigated here, the explicit tau-leaping method is the 
most efficient for this system. 100,000 simulations of the explicit tau-leaping method for the 
small stepsize r = 0.1 took 6.21 seconds with small distance as ones of fully implicit TR-TR 
results for the stepsize r = 0.4. All accelerated methods show efficiency (at least 10 times 
faster) compared to the SSA that took 683 seconds for 100,000 simulations. 

6.3. The ELF System 

We now consider a more complex system containing 8 species and 12 reactions 
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to evaluate the accuracy of the proposed tau-leaping methods. We use the initial conditions 
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and parameter values given in the literature 13(. The chemical reactions, propensity func- 
tions, and initial values are listed in Table 
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Figure 6: The histograms of X$ at the final time obtained with different, fixed stepsizes for the ELF system 
(Table El). Each histo gram uses 100,000 samples. 



We consider the simulation time interval [0, 3] seconds, and perform 100,000 independent 
runs with the Gillespie SSA and with each one of the accelerated methods. The histograms 
of X 5 and X\ concentrations at the final time are presented in Figures [6] and d respectively, 
for different fixed time steps between r = 0.04 and r = 0.005 seconds. Figure [6] shows a 
similar qualitative behavior as in the previous stiff examples. For a large stepsize r = 0.04 
seconds, the histograms produced by the fully implicit BE-BE and BE-TR methods exhibit 



Table 5: List of reactions and propensity functions for the ELF system. 





Reaction 


Propensity 


Rate constant 




Species 


Initial value 


Ri 


E A -> E A + A 


oi = ci [E A ] 


Cl = 


15 




A 


2000 molec. 


R2 


Eb — > Eb + B 


a 2 = c 2 [E B ] 


c 2 = 


15 


x 2 


B 


1500 molec. 


R 3 


E A + B -> E A B 


a 3 = c 3 [E A ][B] 


C3 = 


0.0001 


x 3 


E A 


950 molec. 


Ri 


E A B — > E A + B 


a a = c±[E A B] 


c 4 = 


0.6 


Xi 


E b 


950 molec. 


R 5 


E A B + B -> E A B 2 


a 5 = c 5 [E A B][B] 


C5 = 


0.0001 


x 5 


E A B 


200 molec. 


i?6 


E A B 2 -4 E A B + B 


a 6 = c 6 [E A B 2 ] 


C6 = 


0.6 


x 6 


E A B 2 


50 molec. 


Rr 




a*; = cj[A] 


C7 = 


0.5 


x 7 


E B A 


200 molec. 


R$ 


Eb+A^EbA 


a 8 = c 8 [.Eb]L4] 


C8 = 


0.0001 


x 8 


E B A 2 


50 molec. 


i?9 


E B A —> Eb + A 


a 9 = c 9 [E B A] 


eg = 


0.6 








-Rio 


E B A + A^ E B A 2 


a±0 = ci0[E B A][A] 


ciO = 


= 0.0001 








Rn 


E B A 2 -> E B A + A 


ail = Cl l[E B A 2 ] 


c x l - 


= 0.6 








Rl2 


B -> 


ai2 = ci2[B] 


ci2 = 


= 0.5 
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Figure 7: The histograms of Xi at the final time obtained with different, fixed stepsizes for the ELF system 
(Table E]). Each histo gram uses 100,000 samples. 



a weak damping effect (small sharp peaks), while the histograms given by the implicit order 
two weak Taylor methods with a = 1.0 exhibit a dispersive effect (broader peaks). Figure [7] 
shows a different behavior. For a large stepsize r = 0.04 seconds the BE-BE, the BE-TR, 
and the implicit order 2.0 weak Taylor with a = 1.0 methods show dispersive behavior (broad 
peaks). Therefore the errors in variance for the ELF system have a complex behavior when 
stepsizes are very large. In Figures [6] and [TJ the histograms given by the fully implicit TR- 
TR method and implicit order two weak Taylor method with a = 0.5 are very similar to the 
exact SSA histogram. If the stepsize r is decreased to r = 0.005 seconds, all approximation 
methods show very good accuracy. 

Figures [8] (a) and (b) show the error in distribution (the distance (16T|) between the SSA 
and each of the accelerated methods' histograms) versus simulation stepsize for the ELF 
system. The y-scale in Figure El (b) is much larger than that of Figure M (a) because the 
number of molecules for X\ is much larger than that of X§ (see the Figures |6] and [7]) . The 
results indicate that, similar to the previous examples, the TR-TR and the implicit second 
order weak Taylor method with the a = 0.5 are the most accurate accelerated methods. 

Figure M shows the relationship between accuracy and CPU time for the different stepsizes 
of the ELF system. The accuracy is measured by the distance (161~]) between the accelerated 
method and the SSA histograms for X 5 , as in Figure |8] (a). 100,000 simulation of the SSA 
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Figure 8: The relationship between the error in distribution (the distance ([FT]) between SSA and each of the 
proposed methods' histograms) and the different stepsizes for X 5 and X\ for the ELF system. 
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Figure 9: The relationship between accuracy and CPU time for X§ of the ELF system 

took 178,364 seconds (approximately 50 hours), while 100,000 simulations of the implicit 
order two weak Taylor method with a = 1.0 and (3 = 1.0 for the smallest stepsize r = 0.005 
took 6,216 seconds (3.5% of the SSA time) and provided an accurate solution (distance value 
is only 0.15). For the largest fixed stepsize r = 0.04 seconds, the fully implicit TR-TR and 
the implicit second order weak Taylor method with the a = 0.5 provide high accuracy and 
high efficiency (only 0.4% of the SSA time). 
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7. Conclusions 

This paper develops new implicit tau-leaping-like algorithms for the solution of stochastic 
chemical kinetic systems. The fully implicit tau-leaping methods, "BE-BE", "TR-TR", 
and "BE-TR" , are motivated by the fact that existing implicit tau-leaping algorithms treat 
implicitly only the mean part of the Poisson process. The newly proposed methods treat 
implicitly the variance of the Poisson variables as well. The implicit second order weak 
Taylor tau-leaping methods are motivated by the theory of weakly convergent discretizations 
of stochastic differential equations, and by the fact that Poisson variables with large mean 
are well approximated by normal variables. 

Theoretical stability and consistency analyses are carried out on a standard test problem 
- the reversible isomerization reaction. The fully implicit tau-leaping methods are uncon- 
ditionally stable; the implicit second order weak Taylor tau-leaping methods with a = 1.0 
are conditionally stable, and with a = 0.5 unconditionally stable. The asymptotic means 
of the solutions given by all proposed methods converge to the analytical mean of the test 
problem. The asymptotic variances of the proposed methods, however, converge to different 
values, as it is also the case for traditional tau-leaping methods. 

Numerical experiments are carried out using the decaying-dimerizing system, the bistable 
Schlogl reaction system, and the ELF system to validate the theoretical results. The accuracy 
of the solutions is evaluated by comparing the probability densities obtained with the new 
methods and with Gillespie's SSA. The numerical results verify that the prosed methods 
are accurate, with an efficiency comparable to that of the traditional implicit tau-leaping 
methods. The theoretical analyses and numerical experiments shows that the fully implicit 
TR-TR and the implicit second order weak Taylor tau-leaping methods with a = 0.5 are 
the most accurate methods for large stepsizes. 
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